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Abstract 

We study the dynamics of a colloidal fluid in the full position-momentum phase space. 
The full underlying model consists of the Langevin equations including hydrodynamic in- 
teractions, which strongly influence the non-equilibrium properties of the system. For large 
systems, the number of degrees of freedom prohibit a direct solution of the Langevin equa- 
tions and a reduced model is necessary, e.g. a projection of the dynamics to those of the 
reduced one-body distribution. We derive a generalized dynamical density functional theory 
(DDFT), the computational complexity of which is independent of the number of particles. 
We demonstrate that, in suitable limits, we recover existing DDFTs, which neglect either 
inertia, or hydrodynamic interactions, or both. In particular, in the overdamped limit we 
obtain a DDFT describing only the position distribution, and with a novel deflnition of the 
diffusion tensor. Futhermore, near equilibrium, our DDFT reduces to a Navier-Stokes-like 
equation but with additional non-local terms. We also demonstrate the very good agreement 
between the new DDFT and full stochastic calculations, as well as the large qualitative effects 
of inertia and hydrodynamic interactions. 



1 Introduction 

Since the experimental discovery of the Brownian motion of pollen particles in water in the 19th 
Centuryfli, the study of classical fluids (systems of particles at sufficiently high temperature so 
that quantum effects can be neglected), such as colloidal systems, has been fundamental not only 
to the development of statistical mechanics starting from the work of Einstein^], Langevin[3] 
and Smoluchowskidj, but also to many other fields in physics, chemistry and engineering, e.g. the 
evolution of microscopy over the last century [51 [51 [7], recent advances in biophysical research[5] and 
the rapidly-growing field of microfluidics[^ fTHl lllj . Colloidal systems, in particular, are versatile 
model ones for both theoretical and experimental scrutiny, e.g. many of the forces governing their 
structure and behaviour govern also the structure and behaviour of matter, whilst the sufficiently 
large physical size of colloidal particles makes them accessible experimentally. The large number of 
particles present in real-world systems translates to high-dimensional mathematical models which 
quickly become computationally intractable. 

Statistical mechanics approaches, such as dynamical density functional theory (DDFT) [12], 
allow the dynamics of systems of arbitrarily large numbers of particles to be studied. However, 
existing DDFTs neglect either the momentum of the colloidal particles [T3], or the hydrodynamic 
interactions mediated through the bath[T^, or both[TS|. Yet, inertial effects are negligible only 
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in the high-friction hniit'16], whilst hydrodynamic interactions are long range [TT]: it is therefore 
unclear that existing DDFTs are sufficiently accurate to model colloidal systems of general appli- 
cability. Here we outline a generalized DDFT formalism which includes hydrodynamic interactions 
in the full position-momentum phase space description. In particular, we recover existing DDFT 
formulations fT3 [13 E] as special cases and we further demonstrate the effects and implications 
of including hydrodynamic interactions and inertia. The new DDFT is validated with stochastic 
simulations. 

The careful inclusion of both inertia and hydrodynamic interactions is a large step towards 
accurate and predictive modelling of physically-relevant systems. The numerical experiments 
presented here for model cases of hard spheres in spherically-symmetric external potentials should 
serve as a solid starting point for a wide range of applications. 



2 The Langevin and Fokker-Planck Equations 

We are interested in systems consisting of a large number N of identical, spherically symmetric 
colloidal particles of mass m suspended in a bath of many more, much smaller and much lighter 
particles. Typically, colloidal particles have sizes ranging from Inm-l/^m, occupying the same 
volume as approximately 10^-10^" water molecules. As such, treating the bath particles exactly 
is computationally prohibitive. However, a typical timescale for a colloidal particle to diffuse a 
distance equal to its diameter is Is, whilst the typical time between collisions of water molecules is 
10~^^s [18]. Hence, if we are interested in timescales significantly larger than this typical collision 
time, we may introduce a coarse-grained model and consider only the positions and momenta of 
the colloidal particles, treating the bath in a stochastic or probabilistic manner. 

This approximation leads to the Langevin[5 (or stochastic Newton's) equations for the 3N- 
dimensional colloidal position and momentum vectors r — (ri, . . . , rjv) and p — (pi, . . . , pat) with 
Ti and Pi the position and momentum of the ith particle: 

S^'' VrnM)-r(r)p + A(r)fW. (1) 

Here, V is the potential, which is generally a sum of an external potential, such as gravity, and 
inter-particle forces, such as electrostatic effects. The first equation shows that the change in 
position is given by the velocity, whilst the second shows that the change in momentum is given 
by the force, the first part of which is the negative of the potential gradient. The final two terms 
model the effects of the bath. The motion of the colloidal particles causes flows in the bath, which in 
turn cause forces on the colloidal particles, called hydrodynamic interactions (HI); see Figure [TJa) 
for an illustrative sketch. The momenta and forces are related by the 3iV x 3A^ positive-definite 
friction tensor T, a generalisation of self-friction tensor 7I, where 1 is the 3N x 3A^-dimensional 
identity matrix and 7 is the friction felt by a single isolated particle. Collisions of bath particles 
with colloidal particles are described by the stochastic forces f , given by Gaussian white noise. The 
strength of these collisions is related to T via a generalized fluctuation-dissipation theorem|19j. 
giving A — {mkBTTy^"^ , where T is the temperature and is Boltzmann's constant. Note that 
here and in the following we have assumed that T is constant in space, i.e. that the solvent bath 
is also a heat bath, equilibrating on a much faster timescale than the particle dynamics. 

Equations ([T]) are normally solved numerically by a time-stepping algorithm such as the Euler- 
Maruyama scheme [20]. If HI are neglected by setting F = 7I then A = {mkBT^y^'^l is known 
explicitly and the numerical solution of ([T]) takes 0{N) operations at each timestep. In contrast, 
including HI requires A to be determined from F, e.g. through Cholesky decomposition, which 
then increases the operations required for the numerical solution to 0{N^). Hence, since N is large, 
HI are often neglected to reduce the computational cost, whereas physically they are important 
in all but very dilute systems as they decay only with the inverse of particle separation [ITl [21]. 
Figure[l]demonstrates the effects of including HI in a very simple system, nine identical hard-sphere 
colloidal particles subject only to a constant vertical force —1 with 7 = 1. If HI are neglected, 
the particles accelerate vertically until reaching terminal velocity. Including HI not only increases 
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Figure 1: Dynamics of hydrodynamically-interacting particles, (a) Bath fluid flows (blue) induced 
by motion of spherical colloidal particles (red) restricted to a plane, (b)-(e): Motion of nine 
identical hard-sphere colloidal particles, coloured by symmetry, under a constant vertical force and 
HI at zero temperature (solution of (jlj). (b) Evolution in the centre of mass frames with (left) 
completely symmetric square lattice and (right) slightly perturbed initial conditions, (c) Vertical 
and radial (from centre particle in horizontal plane) positions and velocities of corresponding 
coloured particles for symmetric initial condition. Black curves correspond to no HI. Colloid 
velocities and bath flows in (d) laboratory reference frame and (e) centre of mass frame. Bath flow 
coloured by vertical position. Other quadrants are symmetric. See also Supplementary Movie 1. 
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the speed of sedimentation, but also introduces complex oscillatory dynamics caused by the flows 
in the bath induced by the particle dynamics, which are circulatory when viewed in the centre 
of mass frame; see Figure [l|e). Note that both self-interactions and lubrication effects have 
been neglected. Whilst the former may be measured experimentally [5], they are negligible on the 
timescales described here. However, neglecting lubrication effects leads to a kink in the radial 
velocity (at around non-dimensional time 30) due to the hard sphere effects. 

When N is large, a particular realization of the dynamics of particles, as given by ([!]), is 
irrelevant as it is practically impossible to know the exact initial conditions and, as can also be 
seen from Figure [l|b) , a small perturbation in initial conditions can lead to vastly different long- 
time dynamics. What is of interest are average quantities over a large number of experiments, or a 
large number of realizations given by Q. Mathematically, this corresponds to the joint probability 
distribution f^-^\r,p,t), which is the probability of finding the particles with positions r and 
momenta p at time t. 

Averaging the Langevin equation ([T]) over the initial particle distribution and the noise leads 
to the Kramers equation (the corresponding Fokker-Planck equation in phase space), a 6N- 
dimensional deterministic partial differential equation for the evolution of the distribution function: 

= (r, p, t) + -p • Vr/(^) (r, p, t) 

m 

-Vry(r,i)-Vp/W(r,p,t) (2) 
- Vp • [r(r)(p + mfcsTVp)/(^) (r, p, t)] 

Note that this transformation does not reduce the numerical complexity, as taking M discretization 
points for each spatial/momentum dimension would require (2M)3^ points in total; the only way 
to solve the Kramers equation in high-dimensions is by using Monte Carlo methods, i.e. by 
solving ([I]). However, it is known rigorously f22] that Z*^^^ is a functional of the one-body position 
distribution p{ri,t) — N J dpdr'/*^^^(r, p, i), where dr' denotes integration over all but ri, the 
position of the first particle. Hence, for any number of particles at fixed time, the system is, at 
least in principle, completely described by a function of only a single three-dimensional position 
variable (this is analogous to time-dependent density functional theory in quantum mechanics |i23|). 



3 Dynamical Density Functional Theory Formalism 

The above motivate the derivation of a dynamical density functional theory (DDFT), i.e. an equa- 
tion of the form dtp{ri,t) = — • a(ri,t, [p]), where a is a functional of p. The earliest and 
most basic DDFTs, as suggested by Evans[2¥ and derived by Marconi and Tarazona[l5 , neglect 
both inertia and HI. An improved version that still neglects inertia but includes two-body HI 
was recently given by Rex and L6wen|T3]. Finally, Archer iM; derived a DDFT which includes 
inertia but neglects HI. All these DDFTs are special cases of our general formulation presented 
here, which includes both inertia and HI. Neglecting the HI terms in this new formulation trivially 
returns Archer's DDFT, whilst considering the high-friction limit recovers the configuration-space 
DDFTsfEl [E]; see later. This, along with the comparisons to the full stochastic simulations 
described below, serves as a validation of our generalized DDFT. 

We denote the position and momentum of the jth particle by Vj and pj and let r(r) = 
7[1 + r(r)], where the HI tensor T is decomposed into 3x3 blocks Fij^lG . Physically, r^j 
describes how the momentum of particle j generates a force on particle i. 

To obtain a reduced model for the Kramers equation, we consider its moments with respect to 
momentum and obtain an infinite hierarchy of equations. The challenge is now to find a suitable 
truncation of this hierarchy; this is analogous to deriving the Euler or Navier-Stokes equation 
from the Boltzmann equation[25j. We truncate the hierarchy at the second equation, giving a 
continuity equation 

= dtpiri,t) + Vr, • (p(ri, t)v(ri, t)) (3) 
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and an equation for the time evolution of the current: 



= at(p(ri,t)v(ri,t)) +7p(ri,<)v(ri,t) 

^- / dr'Vr,V^(r,t)p(^)(r,t) (V) 

J 



m 



m 



^ / dpdr'fy(r)p,/W(r,p,t) (H) 

+ Vr,./dp,Pi^/«(r„p„t). (K) 
J 

Here p^^\r,t) = J dp/'^^)(r, p, i) is the iV-body position distribution, Z*^^' (ri, pi, t) = 
J dp'dr'/(^)(r, p, is the one-body phase-space distribution and p(ri, i)v(ri, = 
N/m f dpdr'pi/(^)(r, p, i) is the current. To close the second equation as a functional of p 
and V, it is necessary to deal with terms arising from the many-body part of the potential in (V), 
HI in (H), and 'kinetic pressure' effects in (K). 

For (V), at equilibrium there exists an exact functional identity 

N J dr'Vr,l/(r)pW(r)=p(ri,t)Vr,^-fcBTVr,p(ri), 

where is the Helmholtz free energy functional, 

T[p] = keT J drip(ri)[ln (A3p(ri)) - l] + T,^c[p] 

+ J drip(ri)Fi(ri), 

where A is the Broglie wavelength, and plays no part in the following. In general, J-cxc (the excess 
over ideal gas term) is unknown but has been well-studied at equilibrium and good approximations 
exist, e.g. Rosenfeld's fundamental measure theory[2ni[Il], which is accurate for hard spheres, and 
mean field theory [T2|, which becomes exact for soft interactions at high densities. Motivated by the 
success of such approximations in DFT, we assume that the same equilibrium identity also holds 
out of equilibrium. This has the further advantage of giving the correct equilibrium behaviour. 

Since HI vanish at equilibrium there exists no analogous identity for (H). For ease of exposi- 
tion, we assume that the HI are two-body: Tij{r) = Sij^^e^i'^ii^ij''^^) + (1 ~ %)Z2(ri,rj). We 
also assume that we may write the two-body reduced distribution [5^ as Z*^^^ (ri, r2, pi, P2, ~ 
/(^^(ri, pi, i)/(^)(ri, pi, i)5(ri, r2, [p]) for a known functional g. Again, this has been well-studied 
at equilibrium [24]. To go beyond this two-body approach it is necessary to approximate higher- 
order reduced distributions as functionals of p. 

As (K) is analogous to the kinetic pressure tensor [^7], there is no reason to expect it to 
be a simple functional of p and v. We are motivated by the local-equilibrium approximation 
widely used in statistical mechanicsPT and note that, for a general f^^\ we have /(^^(ri, pi, = 
f^^\ri, pi,t) + /ncq(ri, Pi, i), where /j^,^"* is the local-equilibrium part, the momentum dependence 
of which is given by a local Maxwellian with mass p{r,t), mean mp(r, i)v(r, t) and variance 

mkBTp{r,t). The corresponding mass, mean and variance of the non-equilibrium part /ncq are 
all zero. The local-equilibrium part describes local quantities such as velocity and temperature, 
whilst the non-equilibrium part determines non-local quantities such as heat flux. 

Combining these three approximations and denoting the material derivative Dt = dt + v(ri, t) ■ 
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Vri gives 



Av(ri,i) + ^^Vr,- / dpiHl^/(i)(r,,p,,i) 
p(ri,t) J ^ 



lVr.^^-7v(r.O (4) 
m op 



7 y"dr2p(r2,i)5(ri,r2, [p])Zi(ri, r2)v(ri, 
7 y dr2p(r2,i).g(ri,r2, [p])Z2(ri,r2)v(r2,t), 



which, along with the continuity equation ([3|, and assuming that the term containing /neq may 
be neglected or approximated (e.g. via a Chapman-Enskog expansion) as a functional of p and 
V, gives a DDFT. Neglecting this term can be rigorously motivated both close to local equilib- 
rium, giving a generalised Navier-Stokes equation with non-local terms, and in the overdamped 
limit [16!;, recovering a new Smoluchowski equation with a novel diffusion tensor; see below. Note 
that Q is non-local in both p and v. These non-local terms, absent from previous DDFTs, 
model important correlations in the positions and velocities of different colloidal particles; we 
now discuss their physical interpretation. The term involving Zi combines with 7V to give an 
effective, density-dependent friction coefficient, whilst that involving Z2 couples the velocities in 
a non-local fashion. Surprisingly, this coupling does not require explicit momentum correlations 
in the two-body distribution. These descriptions correspond to the physical interpretation of the 
associated friction tensor submatrices. Neglecting these terms and setting /neq — recovers a 
previous DDFTfM]. Furthermore, setting 7 = gives a DDFT for simple fluids, although in this 
case it is harder to justify neglecting the non-equilibrium term. 



4 Verification and Numerical Experiments 

We now describe three numerical tests of our proposed DDFT. In all cases we solve ([3| and Q for 
hard spheres of diameter 1, and choose m ^ kgT = 1, which is equivalent to non-dimensionalizing 
the equations. We set fil\ = and choose J- to be the hard-sphere FMT functional 26 . Further, 
we choose g to be the simplest possible pair correlation function g{ri, r2, [p]) = 1 for |ri — r2| > 1 
and zero otherwise, which descibes the hard-sphere excluded volume. An extension would be 
to use the analytic correlation function based on the Percus-Yevick equation [55], but the above 
approximation suffices for the systems studied here. 

In Figure |2] we compare the mean radial position and velocity of 50 particles, with friction 
coefficient 7 = 6. The particles begin at equilibrium in a radially-symmetric external potential 
V{r; 3) with minimum at radius 3. At time 0, the potential is instantaneously switched to V{r; 0) 
with minimum at radius 0, which is then instantaneously switched back to V{r] 3) at time 0.5. The 
choice of 50 particles is a compromise between requiring a large number of particles in order to 
overcome the differences between the canonical ensemble stochastic and grand canonical ensemble 
DDFT models, and computational complexity of the stochastic simulations. We show four pairs of 
computations, each containing the solutions of a DDFT (lines) and the corresponding underlying 
stochastic equation (squares) . The first pair (blue) includes both inertia and HI and compares the 
solution of our generalized DDFT ([3| and Q to the Euler-Maruyama[20_ solution of the stochastic 
equations ([T]). The second pair (red) are the same simulations, but when HI are neglected by setting 
r = 7I; this is the special case derived by Archerfl^. As far as we know, these are the first such 
numerical implementations of a phase space DDFT. In both cases, the agreement between the new 
DDFT and stochastic simulations is very good. The quantitative differences for the case including 
HI are due to an enforced difference in choice of friction tensor; see Materials and Methods. In 
addition, the effects of including HI are quite striking. In particular, HI appear to increase the 
effective friction of the system and damp the dynamics. 
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Time 

Figure 2: Mean radial positions and velocities given by (lines) DDFT and (squares) stochastic 
equations. Full phase space with (blue) and without (red) HI from DDFT ^ and Q and stochasics 
([T]). Overdamped limit DDFTJ3J and stochastics|I9j with (green) and without (purple) HI. The 
system contains 50 particles, has 7 = 6, and starts at equilibrium in potential ^ with = 3, 
which is instantaneously switched to tq = at time zero and then back to tq = 3 at time 0.5. See 
Materials and Methods for details. 

The remaining two pairs of simulations are restricted to position space via the high-friction 
approximation, neglecting inertial effects. The DDFTs both with[T3] (green) and without (T5) 
(purple) HI are compared to the Ermak-McCammon|19j solution of the corresponding position- 
space stochastic equations. Whilst the agreement between DDFT and stochastic simulations is 
again very good, it is important to note that neglecting inertia leads to a qualitatively different 
behaviour of the system. In particular, it gives a kink in the mean position, whilst including 
inertia gives smooth curves with a delay before the mean velocity changes sign. Once again, HI 
are seen to significantly damp the dynamics. 

In Figure [3] we show the evolution of the position distribution and velocity for the same system 
of particles, but we now only switch between potentials once, at time 0, and allow the particles 
to evolve towards equilibrium. Once again we compare the solutions of our generalized DDFT, 
both with (blue) and without (red) HI to the corresponding stochastic simulations. Again, the 
agreement between the DDFT and stochastic simulations is very good. The small differences in 
the position distribution near the origin (where the particle density is higher) are likely due to the 
choice of correlation function, which is less accurate at higher densities. Note in particular that 
the inclusion of HI dramatically slows the build-up of particles near the origin. Having verified 
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Figure 3: Radial particle distribution and velocity given by (smooth curves) solution of DDFT 
(H and (|4| versus (noisy curves) stochastic equations ([T]) with (blue) and without (red) HI. Also 
shown is one representative stochastic realization, at times 4, 6, 8 and 10, with particles coloured 
according to their radial position (purple denoting a radial position less than 2, green otherwise). 
The system contains 50 particles, with 7 = 6, and starts at equilibrium in potential ^ with 
To = 5, which is instantaneously switched to rg = at time zero. See Materials and Methods for 
details and also Supplementary Movies 2 and 3. 



our DDFT by comparison to stochastic simulations. Figure |4] depicts the evolution of the position 
distribution and velocity for the DDFT ([s]) and Q with 500 particles, for which the stochastic 
equations are computationally very costly. In this case the effects of including HI are even more 
dramatic, giving qualitatively different behaviour of both the mass distribution and velocity for 
the DDFTs with and without HI. This size-dependence shows that HI must be carefully considered 
in any DDFT used to model macroscopic numbers of particles. 

5 Connection with a Generalized Navier-Stokes Equation 

For the remainder of this discussion we restrict our attention to a potential which is a sum of one- 
and two-body terms and discuss two limits of Q. Close to local equilibrium, we may expand /^^^ 
and /'^^ as Taylor series in [37], obtaining a generalized compressible, non-local Navier- 

Stokes-like integro-differential equation: 

pA.V = r^V^^V + (C + |77)Vr,(Vr, • v) + p G {[p], [v]). 
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Figure 4: Radial particle distribution and velocity given by solution of DDFT ^ and Q 
with (blue) and without (red) HI. The system contains 500 particles, with 7 = 10, initially at 
equilibrium in potential (|6| with vq = 8, which is instantaneously switched to tq = 4 at time 
zero. See Materials and Methods for details and also Supplementary Movies 4. Also shown are 
the number of particles with radii between and 6. 
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where v = v(ri,t), p — p{ri,t) and ^([/o],[v]) is the right hand side of Q. The first three 
terms are standard in the compressible Navier-Stokes equation. However, the viscosities ry and 
C are given by integrals involving the two-body potential and Taylor expansion coefficients of 
/(^) and f'-^\ Hence, the above equation is not amenable to a straightforward numerical solution. 
This restriction is not due to including HI; the same integrals occur when the statistical mechanics 
formulation of the pressure tensor is derived for a simple fluid[27^. The remaining terms (contained 
in Q) are a position-dependent pressure- like term which depends on the gradient of the chemical 
potential /i — 6J-'/Sp, and which is non-local in p, a friction term, and two additional non-local 
terms, the physical interpretation of which was discussed above. 



6 The Overdamped Limit 

As discussed above, most previous DDFTs have been formulated in the high-friction or over- 
damped regime. In such cases, it is assumed that the friction with the bath causes the momentum 
distribution of the colloidal particles to equilibrate on a much shorter timescale than their position 
distribution. Hence one is interested experimentally only in the position distribution whilst the 
momenta, and therefore inertia, may be neglected. In this regime, we have a non-dimensional 
parameter e — ksT / ra'^~^ <^ 1, where L is a 'typical' length scale of the system. De- 
noting a Maxwellian momentum distribution by M(pi) — exp{—\pi\^/{2mkBT))/(2mkBT'K)^^'^, 
we find rigorously p!5] /(^^(ri, pi, t) = M{p){p{ri,t) + ea(ri,t) • pi -I- 0{e^)) for some function 
a; in particular J dpi(pi pi)/'"'^^(ri,pi,t) = 0{e^). For ease of presentation, we assume 
Z2 — (the case Z2 7^ in the overdamped limit is rigorously analysed in reference |16j ) and 
let p(^)(ri, r2, t) — p(ri, i)p(r2, t).g(ri, r2). Then the mass distribution p satisfies a Smoluchowski 
equation 



a,p(ri,0 = Vr, •(D(ri,[p])[^Vr,l/i(ri,<Mri,t) (5) 

+ Vr,p{ri,t) + j^ J dr2p(2)(ri,r2,i)Vr,V2(ri,r2)]), 

where we have a novel, density-dependent diffusion tensor: 

D(ri,[p])-^ 1+ / dr25(ri,r2V(r2,t)Zi(ri,r2) . 

I J J 

Note that since Q is non-linear in p, it is not the Fokker-Planck equation of a Langevin equation. 
Further, ^ can be transformed into a DDFT using standard arguments, but it differs from that 
derived by Rex and Lowen[T3] under analogous assumptions. This is due to the non-commutativity 
of taking the overdamped limit and integrating over all but one particle. 

In contrast to the standard definition, e.g. in reference |T3j, D is a non-local functional of 
the density and thus is also implicitly time-dependent. Surprisingly, this arises even though the 
friction tensor is time-independent. 

The resulting DDFT has the capacity to model significantly more complicated dynamics in 
a much wider range of systems. Previous phenomenological attempts at including a density- 
dependent diffusion coefficient in DDFTs do not correctly take into account the form of the 
diffusion tensor 1 29 . 



7 Conclusion 

Through the inclusion of both inertia and HI, the new DDFT proposed here is general and can 
readily be applied to model a wide spectrum of real-world problems. Some representative appli- 
cations in which we expect the new DDFT not only to provide accurate predictions but also help 
us elucidate the associated underlying phenomena include: (i) the rate, magnitude and pattern 
of aerosol and nanoparticle deposition in the respiratory system, where including inertia quali- 
tatively changes predicted micro- and macroscale deposition, and HI with the walls affect flow 
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patterns; (ii) transport and coagulation of nanoparticles in pulsatile and oscillatory systems, such 
as oscillatory flow mixing, where both inertial and HI effects can quantitatively change the parti- 
cle dynamics (as with the oscillatory system in Figure [2]); and, (iii) cloud formation and aerosol 
production of fine powders, where particle inertia is responsible for preferential concentration and 
HI strongly influence clustering effects. Furthermore, there are many promising extensions to the 
modelling approach proposed here, e.g. to self-propelled particles, modelling bacteria; multiple 
particle species; anisotropic particles; and the inclusion of an external flow field, as would be 
required in modelling blood and drug-laden nanoparticle movement in blood. Similar modelling 
approaches should also be highly relevant in granular media, ion transport, and other multi-phase 
systems. 

8 Methods 

Potential We choose an external potential 

y(r;ro) 0.1(1 - /i)r2 + 3/i - 10 exp[-(r - ro)V4], (6) 

which consists of an asymptotically quadratic confining term, with h a smooth cutoff function 
h{r) = [erf((r 4- ro)/2) — erf((r — ro)/2)]/2, a constant term on the cutoff region, and an inverted 
Gaussian centred at rg. 

Solution of stochastic equations The stochastic equations ([T]) are solved using a standard 
Euler-Maruyama scheme with 10^ time steps. Each simulation is averaged over 5000 runs, the ini- 
tial conditions being chosen via slice sampling the (unnormalized) equilibrium A^-body distribution 
/(^)(r, p) = exp{—V{r)/kBT)exp{—\p\'^/{2mkBT)). The hard sphere potential is approximated 
via a slightly softened, differentiable one 13 . 

Solution of DDFT equations We assume that p(ri,t) and v{ri,t) are functions of only |ri|, 
producing a ID problem. We use spectral methods[3U], appropriately extended to integral op- 
erators, and a fifth order implicit Runge-Kutta method with step size control '3 T. The initial 
condition is obtained by solving the equilibrium DFT equation _24J. 

Friction and diffusion tensors Since T in ([T]) must be positive-definite, we use the inverse 
of the Rotne-Prager approximation with hydrodynamic diameter an — 0.5 [13J. This choice of 
(Th justifies both neglecting lubrication forces and using a two-body expansion. For the DDFT 
(H) and (|4| we use the 11-term two-body expansion given by Jeffrey and OnishiiSl]. Whilst not 
strictly equivalent, the two descriptions are similar when the particles are well-separated in terms 
of the hydrodynamic diameter. In the overdamped limit, the Rotne-Prager approximation to the 
diffusion tensor is used for both stochastic and DDFT calculations. 
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